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Abstract 

This review is a short introduction to numerical hydrodynamics in a cosmological 
context, intended for the non speciahst. The main processes relevant to galaxy formation 
are first presented. The fluid equations are then introduced, and their implementation in 
numerical codes by Eulerian grid based methods and by Smooth Particle Hydrodynamics 
is sketched. As an application, I flnally show some results from an SPH simulation of a 
galaxy cluster. 



1 Gas Physical Processes 



In current cosmological scenarios, the main matter components of the universe are some non 
baryonic Dark Matter (DM), which constitutes most of the universe, and a mixture of primordial 
gas (H, He). The DM component is decoupled from the rest of the universe, and interacts only 
through gravity. The gas component instead can be heated and cooled in several ways, and 
the physics involved is more complicated than in the pure gravitational case. Fortunately, 
while gravity is a long range force, hydrodynamic processes are only important on relatively 
small scales, so that on scales larger than a few Mpc the dynamics of structure formation can 
be studied with good accuracy even neglecting the gas component. Gas dynamics, and the 
related radiative processes, are instead fundamental on smaller scales, e.g. in the formation 
of galaxies, and in linking the matter distribution of the universe to the light distribution we 
actually observe^^ . In what follows I list and discuss very briefly the main gas processes relevant 
to galaxy formation. 

1.1 Heating processes 

Adiabatic compression is the easiest way to heat a gas. By the First Law of Thermody- 
namics, compression work is converted into internal energy: dQ — dU + pdV — 
dU = -pdV. 

Viscous heating is due to the small internal friction (viscosity) present in real gas. Velocity 
gradients in a gas cause an irreversible transfer of momentum from high velocity points to 
small velocity ones, with conversion of bulk velocities into random ones, i.e. into heat, and 
generation of entropy. In the context of galaxy formation viscous heating mostly occurs 
in shocks, which are discontinuities in the macroscopic fluid variables due to supersonic 
flows. These arise for example during gravitational collapse, or during supernova (SN) 
explosions. 

Photoionization takes place when atoms interact with the photons of some background of soft 
X-rays, or UV radiation emitted by QSO or stars: e.g. ^ + H ^ e~ + H^. Observations 
show that at high redshift [z >2) hydrogen in the IGM is indeed ionized (Gunn-Peterson 
effect^)), and although the origin is not clear, this is thought to be caused by some early 
generation of QSO or massive stars. The photoionizing spectrum is usually approxi- 
mated by a power law, with flux J(z/) oc (t//z/£,)"", where vl is the Lyman-a frequency, 
corresponding to the hydrogen ionization energy, 13.6 eV. 

1.2 Cooling processes 

Gas cooling is the key to galaxy formation. In fact, in our current understanding of structure 
formation, the dark matter component of the universe flrst undergoes gravitational collapse, 
forming dark matter halos. These provide the potential wells into which gas can fall and heat 
up by shocks, then immediately cool and form cold, dense, rotationally supported gas disks. In 
these disks conditions are favourable to trigger star formation, and eventually give rise to the 
galaxies we observe today"^-*. The following cooling mechanisms are important in a cosmological 
context. 

Adiabatic expansion is the opposite process of adiabatic compression, with conversion of 
heat into expansion work. 



Compton cooling is electron cooling against the Cosmic Microwave Background (CMB) 
through inverse Compton effect: 7 + e~ —>■ 7 + e~. The condition for this is that 
the temperature of the electrons is higher than the CMB temperature: Tg > T^- The 
net energy transfer depends on the two densities, and on the temperature difference: 
dE/dT oc UeP-yiTg — T^). Since E oc ngTg, the cooling time: tcooi = E/E is in this case 
oc . The CMB photon density decreases like (1 + zY due to the expansion of the 
universe, therefore Compton cooling is only important at high redshifts (typically z >8), 
when the cooling timescale is smaller than the Hubble time. 

Radiative cooling is the most relevant mechanism for the cooling of primordial gas. It is 
caused by inelastic coUisions between free electrons and H, He atoms (or their ions). 
Assuming that the gas is optically thin and in ionization equilibrium, the cooling rate 
per unit volume may be written dE/dt = A{p,T) = rienifiT), where rig and are 
the number densities of free electrons and of atoms (or ions), and f{T) is called cooling 
function. The main processes are: 

• Collisional ionization: the inelastic scattering of a free electron and an atom (or 
ion), which unbinds one electron from the latter, e.g. e~ + if — > + 2e". The net 
cooling for the system is equal to the extraction energy. 

• Collisional excitation + line cooling: the same situation as above, but the atom 
is only excited, and it then decays to the ground state, emitting a photon, e.g. 
e~ -\- H ^ e~ -\- H* e~ + if + 7. This is the dominant cooling process at low 
(10^ K ^10^ K) temperatures. 

• Recombination: e.g. e~ + — > if + 7. A free electron is captured by an ion and 
emits a (continuum) photon. 

• Bremsstrahlung: free-free scattering between a free electron and an ion, e.g. e" + 

if"*" — ^ e~ + if"*" + 7. Its cooling rate grows with the temperature: dE/dt oc T^/^; 
therefore, bremsstrahlung is the dominant cooling process at high (T >10^ K) tem- 
peratures. 

Figure 1 shows the coohng and heating functions in different cases. 
1.3 Other processes 

Besides the heating and cooling mechanisms listed above, other processes may be relevant in a 
galaxy formation scenario. Among them: 

Thermal conduction is direct transfer of heat from regions at high temperature to regions 
at lower temperature, due to the energy transport of diffusing electrons. The induced 
change in internal energy per unit volume is dE/dT — V ^reVrj, with (positive) thermal 
conductivity k — k{T,p). 

Radiation transfer is important in an optically thick medium. Photons are absorbed by gas 
clouds, thermalized by multiple scattering and reemitted as Black Body radiation. This 
process depends on the optical depth Tppt of the gas. 

Star formation follows gas cooling as the next natural step in modeling galaxy formation. 
Our understanding of the detailed physics of star formation is still rather poor, so what 
is usually done is to use some empirical prescription to characterize the gas which is 




Fi gure 1; Cooling and heating functions. Left panel: cooling function A(T)/n'fj versus temperature, for 
a primordial gas. The different contributions to radiative cooling, and the total cooling curve are indicated. 
No photoionizing background is assumed. Right panel: a heating term T/njj is added, due to a photoionizing 
UV background with spectral index a = 1.5. Its effect is both to change the ionization equilibrium (and so to 
change A{T)/n'jj), and to heat the gas. Thin lines correspond to a gas density equal to the mean background 
density at redshift z = 5; thick lines correspond to a density 200 times larger. For low density gas the effect of 
a UV background is dramatic, and line cooling is suppressed. The temperature T^q where | A — r| drops to zero 
is called equilibrium temperature. At T > Teq the gas is cooled, at T < T^q it is heated. This figure was kindly 
prepared by I.Forcada using the atomic rates provided by T.Abel. 

supposed to turn into stars. A good example of recipe^^ is to require: i) a, convergent gas 
flow: V -v < 0; ii) a, Jeans' instability criterion: the free-fall time of the gas cloud be less 
than its sound crossing time; in) a minimum number density of H atoms, e.g. uh > 0.1 
cm~^; iv) minimum gas overdensity, e.g. pg ^py = 178pg, where pv is the virial density 
of the spherical collapse model, and pg is the mean background gas density. If all these 
conditions are satisfied, the gas will form stars at a formation rate similar to the one 
observed in e.g. spiral galaxies. 

Assuming some Initial Mass Function for the stars so formed, it is possible to compute the 
fraction of gas forming massive stars (M > 8Mq). These stars will explode as Type II SN, 
each one releasing 10^^ erg of energy back in the ISM, causing new shocks and gas heating, 
and triggering further star formation. It is also possible to include in this recipe gas release 
and metal enrichment from SN explosions. Unfortunately, at this point the physics of star 
formation is still poorly understood, and the resulting scenarios depend very much on the 
kind of assumptions and modeling made. 



2 Fluid Equations 

For our purpose it is useful to treat the gas as a continuum, and to resort to a hydrodynamic 
description. The fluid equations express conservation of mass (continuity equation), of momen- 
tum (Euler equation) and of energy. We also need an adiabatic state equation: ds/dt — or 



P = P(P) T). These equations may be written as 



|^-,V., (I, 

dv . . , . 

p— = — yp + viscosity terms — pV<P, [2) 

(JjO 

= —pV ■ V + viscosity terms + V ■ f^VT^ + (Q — A) , (3) 

p = p{p,T). (4) 

Equation is given in terms of the specific internal energy e. In Equations (|l|) to (|]) the Ihs 
and the first term on the rhs constitute the usual fiuid equations for a perfect adiabatic gas. The 
extra terms on the rhs are the nonadiabatic terms introduced in the previous Section. Gravity is 
included in the Euler equation via the gravitational potential $. Artificial viscosity terms must 
be included to enable the numerical treatment of shocks and entropy production, processes 
that do not exist in a perfect adiabatic gas. The terms Q and A denote respectively the heat 
sources: photon absorption (e.g. photoionization) and energy feedback from SN, and the heat 
sinks: Compton and radiative cooling, and other photon emission processes. Including all these 
processes (and others, e.g. magnetic fields) in a recipe for galaxy formation is not an easy task. 
Limits arise both from the computational limitations of present day machines, and from our 
ignorance of the physics involved (e.g. the epoch and spectrum of a photoionizing background, 
or the mechanism and role of energy feedback from stars and SN). Moreover, some processes 
naturally go with others, so that if cooling is implemented, star formation and energy feedback 
should also be modeled. For these reasons, different physical processes may or may not be 
taken into account in different computations. For example, current hydrodynamic simulations 
of structure formation sometimes include cooling processes, less often photoionization and star 
formation. On the other hand, thermal conduction and radiative transport have been so far 
neglected. The former may play a role in the central part of a system, but its efficiency 
would depend on the presence of (unknown) magnetic fields. Ignoring the latter is probably a 
good approximation everywhere except in very dense, optically thick regions. Galaxy cluster 
formation is a neat application of these methods, because cooling and star formation are less 
crucial here than in galaxy formation, so they can be ignored to first approximation. An 
example of hydrodynamic simulation of the formation of a galaxy cluster is briefiy presented 
in the last Section of the paper. 

2.1 Eulerian Methods 

There are two basic ways to numerically solve the above set of equations: Eulerian and La- 
grangian. Eulerian methods solve the fiuid equations on a discrete grid fixed in space. In the 
traditional formulation, a Taylor expansion of the terms is used to build a smooth solution 
of the differential equations across different grid cells, as follows. Let / be a one dimensional 
scalar field, and F its fiux; the conservation equation for F is 

df _ OF 
dt dx ' 

if we Taylor expand f{x,t) in time: 

fix, t + dt) = fix, t) + -^dt + ^-g^dt^ + Oidt^) (6) 



and insert Eq. (H) into Eq. (P), we can substitute all time derivatives with spatial ones: 



f{x,t + dt) 



, OF ^ Id 



dF OF 
dx df 



dt^ + 0{dt^ 



(7) 



This is the time evolution equation for /(x, t), which is discretized and solved on a grid, to first 
or second order accuracy. Macroscopic discontinuities in the flow, like those caused by shocks, 
are mimicked by smooth solutions, obtained introducing explicitly, in the fluid equations, the 
artificial viscosity terms mentioned above. The underlying idea is that in a real fluid shocks 
are smooth solutions if seen at small enough scale. 

More recent techniques, named shock capturing schemes or Riemann solvers, use a different 
approach. They incorporate in the method the exact solution of a simple nonlinear problem, the 
Riemann shock tube. This solution describes the nonlinear waves generated by a discontinuous 
jump separating two constant states. The fluid flow is then approximated by a large number 
of constant states for which the Riemann shock tube problem is solved. This automatically 
leads to an accurate approximation of both smooth solutions and of shocks, with no need of 
explicitly introducing viscosity terms in the equations. This category of techniques includes the 
Piecewise Parabolic Method^-* (PPM) and the Total Variation Diminishing^^ scheme (TVD). 



2.2 Lagrangian methods: SPH 

Smooth Particle Hydrodynamics^-'''^'' (SPH) is the most commonly used Lagrangian method in 
dissipative simulations of galaxy formation. By analogy with A^-body codes, where a collision- 
less fluid is represented with a set of discrete particles, SPH also uses particles to describe the 
evolution of a gas fluid. Each particle i is assigned a position r^, a velocity Wj, a density pi 
and a specific internal energy e^, and the fluid equations are solved at the particle's position, 
replacing the true fields with smoothed estimates, which are obtained as local averages of the 
particles' properties. For example, for a scalar field /, the smoothed counterpart (/) is: 

(/(f)) = j d\f{u)W{r- u- h); (8) 

W{r — u] h) is the smoothing kernel, strongly peaked at zero, so that 

lim,^o(/(r-')) = / d\f{u)5D{r- u) = f{r), (9) 

with 5^) a 3 dimensional Dirac delta; h is called the smoothing length. Usually the kernel is 
spherically symmetric: W = W{\f— u|; h), but anisotropic kernels have also been proposed^-*. 

In practice, (/(f)) is evaluated discretely at each particle's position. Defining the particle 
number density at fj as {n{fi)) = p{ri)/mi, in the discrete limit Equation (|^) becomes 

^ m- 

ifi^) = E —m)Wi\f- f,|; h); (10) 
j=i Pj 

for example, for / = p this reduces to 

N 

{p{f)) = Y,mjW{\f-f^\;h). (11) 
i=i 

Typically, in SPH codes the kernel is nonzero only for e.g. |f — u| < 2h, and the smoothing 
length can be varied to keep the summations to the iV = 30 — 50 closest neighbors. 
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Figure 2: Density (left) and temperature (right) radial profiles for an SPH simulation of a galaxy cluster. 
The cluster has a virial mass of 7 x and is resolved by roughly 12000 particles for each species. The 

force resolution is 25 kpc. More details on the figure are given in the main text. 



Following this recipe, one can write the fluid equations: mass conservation is automatically 
satisfied; one form of the Euler and energy equations is, in the adiabatic case''-': 

V,W{\r^-f,\;hy, (12) 
1^ = 3 E ^Av^ - V,) ■ V.W{\n - rSI; h). (13) 

^'^ Pi j=l 

Nonadiabatic terms are introduced in the same way; in particular, shock heating is allowed by 
adding artificial viscosity terms. 

Hybrid methods have also been proposed^°\ that make use of a grid to solve the fiuid 
equations, but are Lagrangian in nature, since the grid itself is deformable and follows the fiuid 
fiow. The quite different approaches of Eulerian and Lagrangian methods, and the variety of 
codes within each approach, make it difficult to compare results coming from different codes, 
and to interpret the comparison. An attempt has however been done^^^ by evolving, from some 
high redshift to the present time, the same initial conditions of a cosmological model, using 
three different Eulerian codes and two different SPH codes; the results seem to indicate that 
Eulerian codes have better resolution on large scales and low p, low T regions, while SPH codes 
can better resolve small scales and high p, high T regions. 



dvj 
dt 



N 

1=1 



2 ' 2 
PI Pi 



3 Applications: Formation of a galaxy cluster 

Some results from an SPH simulation of a galaxy cluster are shown in Figure 2, as an ex- 
ample of application of the theory summarized in the previous Sections. The gas processes 
modeled are adiabatic heating and cooling, and viscous heating, while other cooling processes, 
photoionization and star formation were neglected, since they are not crucial to this particular 
problem. 



The simulation has an Einstein-de Sitter background universe, with Hq = 50 km s~^ Mpc~^, 
and with scale free power spectrum of perturbation P{k) oc k~^. The left panel shows the gas 
density (points) , and the DM and gas mean density profiles (black and grey curve respectively) , 
as a function of the distance r from the cluster center; the vertical bar indicates the virial radius 
of the cluster. Note how the gas is less centrally concentrated than the DM. The density spikes 
in the gas profile are infalling substructure. The right panel shows the gas temperature and the 
mean temperature profile at different radii: the cluster is not isothermal, and its temperature 
drops by a factor of five from the center to the virial radius. During collapse the gas is shock 
heated to about 10^ - 10^ K as seen both in the main system and in the infalling substructures. 
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